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Interesting and recommended websites: 

www.atm.org.uk/mti/2 1 6/data- from- web.html 

A rich resource for data (UK Census, Australian Census and many many others) 

Websites for Dynamic Systems (Prof. Ossimitz, University of Klagenfurt) 
wwwu.uni-klu.ac.at/gossimit/linklist.php (in English) 

wwwu.uni-klu.ac.at/gossimit/sdyn.sdlv.htm Online Kurs zur Systemdynamik 

wwwu.uni-klu.ac.at/gossimit/sw/sdsw.htm Software for Dynamic Systems 

Do you like graphics? Do you want to learn about hyperbolic 
Geometry? Do you like fractals? Then go to: 

www.spektrumverlag.de/sixcms/media.php/767/master2007- 
5-l-small2.pdf 

www-mlO.ma.tum.de/bin/view/MatheVital/IndrasPearls/ 

(Interactive Materials - Cinderella) 

www.pdfone.com/kevword/the-world-of-hyperbolic- 
geometry.html 

Download the pdf-book (and many others) 



Recommended Reading: 

D. Mumford, C. Series, D. Wright: Indra 's Pearls, Cambridge University Press, 2006 


Wonderful Fractals can be found at: 
http://users.telenet.be/ios.leys/ 

Fractals by the famous Fractal's artist Jos Leys (Galleries of Ultra- 
fractal images and animations) 

https://www.fractalus.com/ 

http://www.ultrafractal.com A fractal generating program 
http://infmitezoom.com A fractal gallery 



And finally download a pdf-book for Linear Algebra from 
http://archives.math.utk.edu/tutorials/linear.algebra.krm/ 


Please share your favourite websites and readings with us, Josef. 
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Dear DUG Members, 

Springtime has come to Austria and it is time to publish DNL#77. This newsletter 
is the first one in DUG's twentieth (20 th ) year of existence. I hope that we will 
celibrate this - remarkable for a journal like this - jubilee at TIME 2010 in Malaga. 
We will hold the DUG general assembly in the frame of our conference and I am 
very happy that some of the first members announced attending TIME 2010. Our 
vice president - Barbel Barzel - will be there giving a keynote. Member #1 - Bern- 
hard Kutzler - will give two lectures and there will be many others. We are very 
proud to announce more than 90 lectures and workshops for TIME 2010. Have a 
look on pages 40-42 and browse the accepted abstracts on the conference website. 

Please note the recommended websites. I checked them all and they all worked. In 
my opinion they are really worth visiting. Indra's Pearls ( The Vision of Felix Klein) 
is a great adventure. 

I'd like to inform you that TI-Nspire CAS version 2 is out. Among other changes 
the most remarkable novelty is the availability of colors for plotting and having a 
kind of text and dialogue boxes for programming and within the notes. You can see 
an example on page 32. Updating is an easy process and works fine. Heinrich 
Ludwig's improvement of Runge-Kutta method is a true gem. His contribution in- 
spired me to deal with Runge-Kutta, too - but only with the "classic" method - work- 
ing with Voyage and NspireCAS. 

Wolfgang Alvermann provided the "Octahedron of Horror". This is the translation 
of a headline of an article in the German journal "Der Spiegel" which read there 
"Das Oktaeder des Grauens". The presented problem was part of a (central) end 
examination which appeared to be too difficult for the majority of the students. 
You are invited to making your own opinion and downloading the Spiegel article (URL 
is given on page 39). 

Best regards until summer 



As announced earlier we will put all proceedings of ACDCA- and DERIVE Con- 
ferences on the RFDZ-website http ://rf dz . ph-noe . ac . at/index. php?id= 1 33 . We 
are starting with TIME 2008. 


Download all Z)/VL-Derive- and Tl-files from 

http : / / www . aus tromath . at/ dug/ 


Attend TIME2010 in Malaga, Spain 

http: //www. time2010 .uma.es 
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Contributions: 

Please send all contributions to the Editor. 
Non-English speakers are encouraged to 
write their contributions in English to rein- 
force the international touch of the DNL. It 
must be said, though, that non-English 
articles will be warmly welcomed nonethe- 
less. Your contributions will be edited but 
not assessed. By submitting articles the 
author gives his consent for reprinting it in 
the DNL. The more contributions you will 
send, the more lively and richer in contents 
the DERIVE & CAS-77 Newsletter will be. 


Next issue: June 2010 

Deadline 1 5 May 2010 

Preview: Contributions waiting to be published 

Some simulations of Random Experiments, J. Bohm, AUT, Lorenz Kopp, GER 

Wonderful World of Pedal Curves, J. Bohm 

Tools for 3D-Problems, P. Liike-Rosendahl, GER 

Financial Mathematics 4, M. R. Phillips 

Hill-Encription, J. Bohm 

Simulating a Graphing Calculator in DERIVE, J. Bohm 
Henon & Co, J. Bohm 

Do you know this? Cabri & CAS on PC and Handheld, W. Wegscheider, AUT 

An Interesting Problem with a Triangle, Steiner Point, P. Liike-Rosendahl, GER 

Overcoming Branch & Bound by Simulation, J. Bohm, AUT 

Diophantine Polynomials, D. E. McDougall, Canada 

Graphics World, Currency Change, P. Charland, CAN 

Cubics, Quartics - interesting features, T. Koller & J. Bohm 

Logos of Companies as an Inspiration for Math Teaching 

Exciting Surfaces in the FAZ / Pierre Charland 's Graphics Gallery 

BooleanPlots.mth, P. Schofield, UK 

Old traditional examples for a CAS - what's new? J. Bohm, AUT 

Truth Tables on the Tl, M. R. Phillips 

Advanced Regression Routines for the TIs, M. R. Phillips 

Where oh Where is IT? (GPS with CAS), C. & P. Leinbach, USA 

Embroidery Patterns, H. Ludwig, GER 

Mandelbrot and Newton with DERIVE, Roman Hasek, CZ 

Snail-shells, Piotr Trebisz, GER 

A Conics-Explorer, J. Bohm, AUT 

Practise Working with times 

Huffman-Code with DERIVE and TI-CAS, J. Bohm, AUT 
Tutorials for the NSpireCAS, G. Herweyers, BEL 
Some Projects with Students, R. Schroder, GER 
Dirac Algebra, Clifford Algebra, D. R. Lunsford, USA 

and others 

Impressum: 

Medieninhaber: DERIVE User Group, A-3042 Wiirmla, D'Lust 1, AUSTRIA 
Richtung: Fachzeitschrift 
Herausgeber: Mag. Josef Bohm 


The DERIVE-NEWSLETTER is the Bulle- 
tin of the DERIVE & CAS- 77 User Group. 
It is published at least four times a year 
with a contents of 40 pages minimum. The 
goals of the DNL are to enable the ex- 
change of experiences made with DERIVE, 
77-CAS and other CAS as well to create a 
group to discuss the possibilities of new 
methodical and didactical manners in 
teaching mathematics. 

Editor: Mag. Josef Bdhm 
D'Lust 1, A-3042 Wiirmla 
Austria 

Phone: ++43(0)6604070480 

e-mail: nojo.boehm@pgv.at 
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Danny Ross Lunsford antimatter330yahoo . com 

Hello Prof. Boehm, 

Are you going to make PDFs of the newsletters in the range 19-current? These are a valuable histori- 
cal resource for serious DERIVE users. 

Even though TI has stopped selling it, I use it constantly. 

-ross in atlanta 

DNL: 

Dear Derivian, 

many thanks for your kind and encouraging words. Yes, I intend to proceed revising and 
converting to pdf-files the DNLs which have not been published on the web so far. The next 
one is ready and I hope that it will appear on our website within the next week. 

I also will update regularly an index for the revised issues and the other ones as well. 

The index can be downloaded from the DUG-site. 

Best regards 
Josef 

Danny Ross Lunsford antimatter330yahoo . com 

Sehr gut :) 

Recently I was working with my old favorite APL2 - 1 have a half-hearted project to implement APL2 
in Derive - some things are done - but I came up against the notion of an empty array with structure, 
e.g. 3 0 2 reshape empty creates an empty array with dimensions 3 0 2. Although Derive's array engine 
is fully capable of handling heterogeneous arrays of any dimension and includes the notion of an 
empty array [], I could not find a natural way to implement structured empty arrays (having one di- 
mension = 0). 

In the process of this, it struck my how much more felicitous is Derive than APL2 in actual use. e.g. 
matrix multiplication in APL2 is very ungainly and always requires an explicit operator +.x to be 
placed between factors. I work in the Dirac algebra of spacetime and with higher Clifford algebras, 
and this seemingly minor issue becomes a huge drain on one's energy because it is almost impossible 
to enter expressions of any complexity without syntax errors. In Derive one proceeds in a way almost 
like on paper! I was thinking of publishing my Dirac algebra work in the newsletter, and later a more 
general article on Clifford algebras abstractly. 

-ross in atlanta 
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Phil Ignacio pdyl971@gmail . com 

Is there a Math wizard that help with this please? 

Thanks 

Phil 


This is my problem: 

Given Y\ and its first derivative as t 2 + 5t and 2^ + 5 and Y 2 and its first derivative t 2 - 5t and 2t- 5. 


The Wronskian is 



^2 


Hk 2 . 


Now assume we <DO NOT> know 7 2? however we <DO> have the values of W and Y x . 

How do I get back to the original Y 2 and its first derivative? I tried with Y 2 = a and Y 2 = b and failed. 


Ignacio Larrosa Canestro ilarrosa0MUNDO-R . COM 

You are looking for a function y(t). You must substitute a = y(t) and b = y'(t). Then, you can 
get a differential equation. You can solve it with DSOLVE1_GEN, by example, see the at- 
tached file 

Ignacio showed a unique solution. I believe that the answer is not unique. However I followed Igan- 
cio's advice, Josef. 


#21: y(x) := 


#22: DET 


t + 5 • t y(t) 

2 • t + 5 y 1 (t) 


10. t 


2 2 
#23: (t + 5 • t) *y r (t) - (2-t + 5).y(t) = 10-t 

2 2 

#24: DS0LVE1_GEN(- (2-t + 5)-z - 10-t , t + 5-t, t, z) 


z + 10. t 


#25: 


= c 


t-(t + 5) 
#26: SOLVE 


z + 10.t 

= c, z = (z = t * (c • ( t + 5) - 10)) 

< t-(t +5) J 

For applying DSOLVE I replace y{t) by z (as a function of t). 

A family of solution functions: 

#27: VECTOR(y(t) = t - (c . (t + 5) - 10) , c p -5, 5) 

#28: [y(t) = - 5 • t • (t + 7), y(t) = - 2-t.(2-t + 15), y(t) = - t-(3-t + 25), y(t) 
2 -t • (t + 10), y(t) = - t-(t + 15), y(t) = - 10-t, y(t) = t-(t - 5), y(t) : 


2 -t , y(t) z t • (3 -t + 5), y(t) = 2.t-(2-t + 5), y(t) 


z 5 • t • (t + 3)] 


2 

As you can see one of the solution functions is y = t - 5t (c = 1). 
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This is in Meta-language (but not yet in DERIVE syntax): 

IF MOD(/c,/)=0 
THEN 
IF vl/c=1 

THEN door will be closed 
ELSE: door will be opened 
ELSE: door remains unchanged. 

In DERIVE syntax it reads: lF(MOD(k,i)=0,lF(vlk= 1,0,1) ,vlk). 

We will have our DERIVE procedure as general as possible. The door numbers k shall run 
from 1 through n (the Dodon case will be n = 100). 

What happens on day / can be described in another vector with vector v as its base: 

VECT0R(IF(M0D(k,i) = 0, IF(v4k = 1, 0, 1), vlk) , k, 1, n) . 

We can check this. Take the state of day #4 (the first 20 doors = dim(v_)) from above, name 
it v_ then our procedure should return the state of the 20 doors on day / = 5: 

v_ := [1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1] 

VECTOR (I F(M0D(k, 5) = 0, IF(v_ =1, 0, 1), v_ ) , k, 1, DIH(v_)) 

k k 

[1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0] 

It makes sense to show the number of the day together with the state of the doors, so we 
include this number and define a function of three variables (/ = number of the day, v = last 
state of the doors, n = total number of the cell doors). 

F(i , v, n) := f-i, VECTOR (I F (M0D(k , i) = 0, IF(v = 1, 0, 1), v ), k, 1, n)l 

L k k J 

As we want to work recursive (the state of one day is always the consequence of applying 
the same rule on the state one day before - state(v,) = rule(state(v/- 1 )). Each recursive proce- 
dure needs an initial element which is in our case obviously the state of the doors on day #1 
- all doors are open = all elements of the vector are 1 : 

start(n) := VECT0R(1, k_, n) 

The ITERATES command is an excellent tool for programming recursive procedures. The 
recursion which creates the vector v, +1 from the vector v, (and the day# /'+ 1 from the day# /', of 
course). 

prison(n, m) := ITERATES(F(n + 1, v, n) , [i, v], [1, start(n)], m - 1) 
prison(20, 10) 

This command should show the state of the first n = 20 doors for the first m = 10 days: 






D-N-L#77 


Roland Schroder: The Birthday of the King 


1 

[ 1 , 

1 , 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1] 

2 

[ 1 , 

0, 

1, 

0, 

1, 

o, 

1, 

o, 

1, 

o, 

1, 

o, 

1, 

o, 

1, 

o, 

1, 

o, 

1, 

0] 


[ 1 , 

0, 

0, 

0, 

1, 

1, 

1, 

0, 

o, 

0, 

1, 

1, 

1, 

0, 

0, 

0, 

1, 

1, 

1, 

0] 

4 

[ 1 , 

0, 

0, 

1, 

1, 

1, 

1, 

1, 

o, 

0, 

1, 

0, 

1, 

0, 

0, 

1, 

1, 

1, 

1, 

1] 

5 

[ 1 , 

0, 

0, 

1, 

0, 

1, 

1, 

1, 

o, 

1, 

1, 

0, 

1, 

0, 

1, 

1, 

1, 

1, 

1, 

0] 

6 

[ 1 , 

0, 

0, 

1, 

0, 

0, 

1, 

1, 

o, 

1, 

1, 

1, 

1, 

0, 

1, 

1, 

1, 

0, 

1, 

0] 

7 

[i, 

0, 

0, 

1, 

0, 

0, 

o, 

1, 

o, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

0, 

1, 

0] 

S 

[ 1 , 

0, 

0, 

1, 

0, 

0, 

o, 

0, 

o, 

1, 

1, 

1, 

1, 

1, 

1, 

0, 

1, 

0, 

1, 

0] 

9 

[ 1 , 

0, 

o, 

1, 

o, 

o, 

o, 

o, 

1, 

1, 

1, 

1, 

1, 

1, 

1, 

o, 

1, 

1, 

1, 

0] 

10 

[ 1 , 

0, 

o, 

1, 

o, 

o, 

o, 

o, 

1, 

o, 

1, 

1, 

1, 

1, 

1, 

o, 

1, 

1, 

1, 

1] 


You may compare with prison(10,20). 

Play a little bit with this function until you get an idea which doors will stay open on Dodon's 
great day. Can you find the rule behind? 

It would be great if we could add a visualisation of the situation(s). Day #5 for the first 30 
doors and the history for the first 40 days could look like this: 

#15: doors (30, 5 ) 



#18: time_doors(40 , 30) 
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We leave this for the reader. The procedures can be found in the file. 
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RK6-A Modern Runge-Kutta-Method 

Heinrich Ludwig, Raubling, Germany 

1. Ein Kurzportrait der Verfahren 

Runge-Kutta-Verfahren sind eine Familie von Losungsverfahren fur gewohnliche Differentialglei- 
chungen (DG1) erster Ordnung / = /{x,y) mit gegebenen Anfangsbedingungen. Sie werden eingesetzt, 
wenn die analytische Losung der DG1 nicht bekannt ist und eine Naherungslosung als Ersatz geniigt. 

Computeralgebrasy steme bringen numerische Verfahren zum Losen von DGln mit, in vielen Fallen 
ein explizites Runge-Kutta-Verfahren (RK). Auch DERIVE stellt mit der Funktion RK ein solches 
Verfahren zur Verfiigung. Es ist das „klassische“ RK-Verfahren, das Runge schon 1895 veroffentlicht 
hat. Heute ist es nur noch von historischer Bedeutung; seine Weiterentwicklungen sind wesentlich 
efflzienter und haben zudem erganzende mitzliche Eigenschaften. Mein Ziel war es, DERIVE mit 
einem modernen RK-Verfahren auszustatten, um fur Simulationen in meinem Physikunterricht ein 
leistungsfahiges Werkzeug zur Verfiigung zu haben. 

Den Algorithmus eines RK-Verfahrens zu erklaren wiirde den Rahmen dieses Beitrags sprengen. Viele 
Bucher, die in die Grundlagen der numerischen Mathematik einfuhren, behandeln das Thema. Beson- 
ders ausfuhrlich ist das Buch von E. Hairer, S. P. Norsett, G. Wanner : Solving Ordinary Differential 
Equations , 2. Auflage, Springer- Verlag, 1993. Ein guten Uberblick fiber die jungere Entwicklung auf 
dem Gebiet bietet der Aufsatz von W. H. Enright, D. J. Higham, B. Owren, Ph. W. Sharp: A Survey of 
the Explicit Runge-Kutta Method , TR 291/94, Department of Computer Science, University of Toron- 
to, 1994 (zur Zeit verfiigbar fiber www.cs.toronto.edu/pub/reports/na/rk.survey.96.ps. gz) [1] . 

Ich habe das Verfahren von Papacostas, Papageorgiou und Tsitouras (SIAM J. Numerical Analysis , 33 
(1996) S. 917-936) fur eine Umsetzung in DERIVE-Code gewahlt. Das Verfahren gehort zu den effi- 
zientesten, die in den letzten Jahren veroffentlicht wurden. Der Algorithmus selbst ist in Ausdruck #7 
formuliert. Er benotigt eine Reihe von Koeffizienten, die in der Matrix A 5 und in den Vektoren B5, C5 
und D5 zusammengefasst werden. 

1. A short description of the methods 

Runge-Kutta-Methods are a family of solution methods for ordinary first order differential 
equations with given initial conditions. They are applied if the analytical solutions of the ODE 
is unknown and an approximative solution is sufficient. 

Computer algebra systems offer numerical methods for solving differential equations which in 
most cases include an explicit Runge-Kutta-procedure (RK). So does DERIVE with function 
RK making such a procedure available. This is the “classic” method published by Runge in 
1885. Nowadays it is only of historic importance, its further developments are much more 
efficient and contain additional useful features. It was my aim to implement a modern RK- 
procedure in DERIVE in order to have a modern and efficient tool for simulations in my phys- 
ics teaching available. 

Explaining the algorithm of a RK-method would break the contents of this contribution. There 
are many books, treating the basics of numerical mathematics dealing with this subject. Es- 
pecially comprehensive is the book written by E. Hairer, S. P. Norsett, G. Wanner: Solving 
Ordinary Differential Equations, 2. Edition, Springer-Verlag, 1993. A good survey over the 
latest development on this field is given by a paper published by W. H. Enright, D. J. Higham, 
B. Owren, Ph. W. Sharp: A Survey of the Explicit Runge-Kutta Method , TR 291/94, Depart- 
ment of Computer Science, University of Toronto, 1994 (available for download at 
www.cs.toronto.edu/pub/reports/na/rk.survev.96.ps.ciz ). [1] 


[1] Sie konnen die pdf-Datei unter downloadbaren Dateien zu diesem DNL finden. 

[1] You can find the pdf-file among the downloadable files connected with this DNL. 
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Value was set on the demand that function RK6 can be used in the same way as function 
RK( r , v , vO , h , n) which is implemented in DERIVE. So it is easy to improve existing dfw- 
and mth-files by only replacing RK by RK6 or by replacing RK by RK6 in the utility file ODEAp- 
proximation .mth and then rename RK6 as RK. Of course, A, B5, C5 and D5 must also be 
included in the utility file. 

Das eigentliche RK-Verfahren besteht nur aus der Schleife innerhalb der Funktion RK6. Das Ergebnis 
eines Schleifendurchgangs dient als Startwert fur den nachsten Durchgang. Nach n Durchgangen gibt 
RK6 sein Ergebnis zuriick. 

Ein Mangel von RK, der notwendigerweise auch RK6 anhaftet, ist, dass die Schrittweite aller RK- 
Einzelschritte dieselbe ist. Die Folge sind von Schritt zu Schritt stark schwankende lokale Diskretisie- 
rungsfehler liber das Integrationsintervall [ vOj.1 , vOtl+h • n], Modeme RK-Module haben die Fa- 
higkeit, die Schrittweite selbsttagig so zu steuem, dass der lokale Diskretisierungsfehler annahemd 
konstant bleibt. Die Verfahren werden dadurch nicht nur erheblich effizienter; bei manchen Anfangs- 
wertproblemen ist die fortwahrende Anpassung der Schrittweite unbedingt notwendig. Aus diesem 
Grund habe ich RK6 um eine Schrittweitensteuerung (step size control) zur Funktion RK6A erganzt. 
RK6A wird an der Stelle erklart, wo sie fur eines der folgenden Beispiele verwendet wird. 

The essence of the RK-procedure is given by the loop within function RK6. The result of one 
run of the loop serves as initial value for the next run. RK6 returns its result after having per- 
formed n loops. 

There is one deficiency of RK and therefore for RK6 as well: the step width is the same for all 
single RK-steps. The consequence of this deficiency are strong step by step changing local 
discretization errors over the interval of integration [vOj.l,vOj.l+h*n]. Modern RK-modules 
have the ability to control the step width automatically in such a way that the local discretiza- 
tion error approximately remains constant. These procedures do not become only more effi- 
cient but for some initial value problems perpetual adaptation of the step width is an absolute 
need. This was reason enough for me to improve RK6 to function RK6A by including a step 
size control. RK6A It will be explained when it is used for one of the following examples. 

2. Der Genauigkeitsgewinn mit RK6 

Die erste Aufgabe fur RK6 sei das Keplerproblem, die Bewegung eines Massepunkts in einem radialen 
l/r 2 -Kraftfeld bei gegebenen Anfangsbedingungen. Fur den Ortsvektor r = [x, y ] des Massepunkts gilt 
die DG1 

r"= f(r) 

wobeiXt:) das Kraftfeld ist, in dem sich der Massepunkt bewegt. In diesem Fall ist das Kraftfeld zeit- 
und geschwindigkeitsunabhangig. Gesucht ist die Bahn ijl). RK-Verfahren losen Differentialglei- 
chungen 1. Ordnung. Deshalb muss die o.g. DG1 2. Ordnung in ein System von gekoppelten DGln 
1 . Ordnung umgewandelt werden. Dabei ist v = [vx, v;v| der Geschwindigkeitsvektor des Massepunkts. 

(1) r' = v 

( 2 ) Y'=J(r) 

Beide DGln lost das RK-Verfahren gleichzeitig, wenn man v und f(r) zu einem 4-dimensionalen Vek- 
tor zusammenfasst. Die ersten beiden Parameter der Funktion RK6 im Ausdruck #8 zeigen, wie das 
gemeint ist. Die Spaltenselektion a [ 2 , 3 ] am Ende des Ausdrucks wahlt aus den Ergebnisvektoren 
der Funktion RK6 die beiden Ortskoordinaten [x, y] = r und lasst die beiden Geschwindigkeitskoordi- 
naten weg. 
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2. The win of accuracy working with RK6 

The first task applying RK6 is Kepler's Problem, i.e. the movement of a mass point in a radial 
Mi 2 field of force with given initial conditions. To the position vector r= [x, y] of the mass 
point applies the following DE 

r"= f(r) 

with f(r) being the force of field where the mass point is moving. In this case the field is inde- 
pendent of time and velocity. We want to find the orbit r[t). RK-methods solve 1 st order differ- 
ential equations. So we need to transform the DE from above which is of 2 nd order into a sys- 
tem of coupled 1 st order DEs. v = [vx,vy] is the velocity vector of the mass point. 

(1) r'=v 

(2) v' = f(r) 

The RK-procedure solves both DEs simultaneously if one combines v and f( r) to one vector 
of 4 dimensions. The first two parameters of function RK6 in expression #8 explain how this is 
meant. Selection of columns 1 1 [ 2 , 3 ] at the end of the expression selects the two position 
coordinates [x,y] = rfor ouput (and for plotting) from the resulting vectors and omits the ve- 
locity coordinates. 


#8: 



r 


bahnl :i 

RK6 



k. 

V _ 


1 

x, 

2 2 3 

7(x + y ) 


ji 

[0, 0.2, 0, 0, 3], , 

50 


50 


U[2, 3] 


1 

y . 

2 2 3 

V(x + y ) 


[t, x, y, vx, vy], 


Dieses spezielle Anfangswertproblem wurde gewahlt, weil es fur t = n eine exakte Losung gibt. Mit 
ihr lasst sich die vom RK-Verfahren gefundene Losung vergleichen und die Genauigkeit des RK- 
Verfahrens bewerten. Wenn die Kraft auf den Massepunkt zum Zentrum [0, 0] weist und seine An- 
fangsgeschwindigkeit nicht zu groB ist, stellt \x(t),y(t)] eine Ellipse dar. Die Anfangswerte 
x(0) = 0.2, j(0) = 0, v x (0) = 0 und v^O) = 3 wurden so gewahlt, dass die Bahnellipse die groBe Halb- 
achse a = 1 hat und die Bahn in der Periapsis beginnt. (Periapsis ist dabei der Punkt auf der Umlauf- 
bahn mit der kleinsten Entfernung zum Hauptkorper und Apoapsis der mit der groBten.) Das Produkt 
aus Schrittweite h und dem Wiederholungszahler n wurde gleich h-n = n gesetzt. Das entspricht der 
halben Umlaufdauer auf einer Ellipse mit Halbachse a= 1. Wenn RK6 gut rechnet, sollte der letzte der 
ermittelten Losungspunkte die Apoapsis der Bahn sein, die um 2 a von der Periapsis entfemt liegt , 
also bei [-1.8, 0]. 

Wie gut die von RK6 ermittelte Losung die Apoapsis trifft, zeigt die folgende Abbildung. Der Fehler 
ist im gewahlten MaBstab nicht zu erkennen. (Er betragt 0,00047.) Zum unmittelbaren Vergleich wur- 
de mit Ausdruck #9 das gleiche Anfangswertproblem auch mit dem „klassischen“ RK-Verfahren ge- 
lost. Der Unterschied der Ergebnisse ist deutlich. 

Wichtig: Die Funktionen RK6 und RK6A miissen durch Approximation vereinfacht werden (Tasten- 
kombination Strg+G). Andernfalls steigt die Rechenzeit erheblich an. Eine exakte Vereinfachung 
bringt schon deshalb keinen Gewinn, weil ein RK-Verfahren ein Anfangswertproblem grundsatzlich 
nur naherungsweise lost. Die Plots konnen auch ohne vorherige Berechnung im Algebrafenster direkt 
gezeichnet werden - wie zB den Ausdruck #9 auf der nachsten Seite. Es muss aber die Option Ap- 
proximate before Plotting 64 aktiviert werden. 
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This special initial value problem was chosen because of its exact result for t = tt. So it is 
possible to compare the approximative solution given by the RK-method with the exact 
solution and we can evaluate the accuracy of this procedure. If the force acting on the mass 
point is directed towards the centre [0, 0] and its initial velocity is not too high then [x(t),y(t)] 
appears as an ellipse. The initial values x(0) = 0.2, y(0) = 0, v x (0) = 0 and i/ y (0) = 3 were cho- 
sen such that the major axis of the orbit ellipse a = 1 and the orbit starts at the periapsis. 
(Periapsis the point on the orbit with minimum distance from the mass M, and aoapsis is the 
point with maximum distance.) The product of step size h and the repetition counter n was 
set h-n = tt. This is equivalent to the half period of revolution on an ellipse with half major 
axis a = 1 . If RK6 does really a good job then the last of the calculated solution points should 
be the aoapsis of the orbit which is at [-1 .8, 0] - in a distance of 2 a from the periapsis. 

The quality of the solution calculated using RK6 meets the aoapsis shows the following fig- 
ure. The error cannot be recognized in the chosen scale. (It is 0,00047.) To have an immedi- 
ate comparison the same IVP was solved by the “classic“-RK-method (expression #9 which 
is part of the utility file). The difference of the results is clear. 


Important: Functions RK6 and RK6A must be simplified by approximation (key combination 
Ctrl+G). Otherwise calculation time would increase enormously. There is no benefit in exact 
Simplification because a RK-method basically solves an IVP only approximatively. It is not 
necessary to do the calculation in the Algebra Window. You can immediately plot the high- 
lighted expression (like #9 below) in the 2D Plot Window but you have to activate there op- 
tion “Approximate before Plotting”. 



f 

r 

1 

1 

bahn2 := 

RK 


vx, vy, 

■ x, -y 

#9: 



2 2 3 

2 2 3 
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Ax + y ) Ax + y ) 
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[0, 0.2, o, o, 3], , 

50 
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1 1 [2 , 


3] 



Wie man sieht, iibertrifft RK6 RK an Genauigkeit bei weitem. Freilich kann man die Schrittweite h so 
weit verringem, dass auch RK keinen so groBen Fehler mehr erzeugt. Eine Verringerung von h kommt 
aber RK6 noch viel mehr zugute als RK, da der Fehler von RK6 ungefahr proportional zu h" abnimmt, 
der von RK nur proportional zu h 4 . Gerade, wenn hohe Genauigkeit erforderlich ist, sind Verfahren 
hoher Konvergenzordnung effizienter. 
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As you can see RK6 is much more accurate than RK. Though one could decrease step size h 
in such a way that even RK does not produce such errors. But then RK6 benefits much more 
of the decreased step size than RK because the error of RK6 declines proportional to h 6 and 
the error of RK declines only proportional to h 4 . Especially when high accuracy is needed 
procedures of high order of convergence are more efficient. 


Bewegung in drei Dimensionen 


Die Physik kennt interessantere Anwendungen fur RK-Verfahren als das Keplerproblem. Wie sieht 
etwa die Bahn eines geladenen Teilchens aus, das sich einem magnetischen Dipol nahert? Ein Beispiel 
dafiir sind die Partikel des Sonnenwindes, die sich unter dem Einfluss des Erdmagnetfelds der Erde 
nahem. Die Bahn der Teilchen ist nicht eben, ihre Berechnung muss in drei Dimensionen erfolgen. 
Ausdruck #10 liefert die rechte Seite der Bewegungsdifferentialgleichung. Man erkennt darin das 
Kreuzprodukt, das fur die Lorentz-Kraft kennzeichnend ist. Im Ausdruck #12 wird die Bahn berech- 
net. 


3. Motion in three dimensions 

In physics we can find applications of RK-procedures which are much more interesting than 
Kepler's Problem. Let’s ask for the orbit of a charged particle approaching a magnetic di- 
pole? One example for this are the particles of the solar wind approaching earth influenced 
by its magnetic field. The orbits of the particles are not plane, their calculation needs three 
dimensions. Expression #10 delivers in its right side the differential equation of the motion. 
One can recognize the cross product which is significant for the Lorentz force. Expression 
#12 calculates the orbit of the particle. 


magnetdipol (w, a, e_, f_, g_) := 

Prog 

e_ := [3-W|2-W|4, 3-wi3-wi4, 2 -wj, 4 a 2 - w;2 A 2 - W|3 A 2] 
f_ := [wi5 , wi6, wi7] 

g_ := CROSS (f_, e_) • a/7(wi2 A 2 + wi3 A 2 + wi4 A 2) A 5 
RETURN APPEND(f_, g_) 
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#12: bahn3 := (RK6(magnetdipol (q , -1), q, [0, 0.1, 0, 1.5, 0, 0, -0.6], 0.01, 


460)) ii [2 , 3, 4] 
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Ein 3D-Graphikfenster zeigt diese Bahn (linkes Bild). Aufschlussreich ist auch das zugehorige Ge- 
schwindigkeitsdiagramm (rechtes Bild). Dazu ersetzt man in #12 die Spaltenselektion 1 1 [ 2 , 3 , 4] 
durch i i [ 5 , 6 , 7 ] . Man sollte fur dieses Bild, die Schrittweite h des RK-Verfahrens verringern und 
die Schrittzahl n entsprechend erhohen, damit die Interpolationssehnen des Graphen keine deutlichen 
Ecken bilden. Der Verlauf des Graphen auf einer Kugeloberflache ist kein Zufall: In einem zeitlich 
konstanten Magnetfeld bleibt die kinetische Energie eines Teilchens und damit v 2 = v x 2 + v y 2 + v z 2 kon- 
stant. 

A 3D-Plot Window shows this orbit (left figure). The respective velocity diagram (right figure) 
is very informative. Replace in #12 selection of columns u[2 , 3 ,4] byii[5,6,7]. One 
should decrease step size h and increase step number n for obtaining a smooth graph. The 
run of the orbit on the surface of a sphere is not pure chance: in a temporal constant mag- 
netic field the kinetic energy of a particle which is v 2 = v x 2 + v y 2 + v 2 remains constant. 



Dieselbe Differentialgleichung mit anderen Anfangswerten liefert eine vollig andere Bahn. Ein Bei- 
spiel dafiir ist die chaotische Bahn eines Teilchens, das im Van-Allen-Gurtel gefangen ist. Ausdruck 
#13 berechnet sie. Wegen seiner reizvollen Erscheinung wird auch der Verlauf des Geschwindigkeits- 
vektors abgebildet (Ausdruck #14 - nachste Seite, rechtes Bild). Aber auch die Bahn selber ist se- 
henswert (nachste Seite, linkes Bild). 

Beide Ausdrucke verwenden vorteilhaft RK6A anstelle von RK6. Die Parameterliste von RK6A unter- 
scheidet sich von der von RK6 bzw. RK. Erlauterungen dazu finden Sie im Anhang. 

The same DE with other initial values leads to a very different orbit. One example for this is 
the chaotic orbit of a particle which is caught in the Van-Allen-Belt. It is calculated in expres- 
sion #13. In addition we show the run of the velocity vector because of its pretty appearance 
(next page, right figure). But the orbit itself is also worth seeing (next page, left figure). 

Both calculations use the advantage of RK6A compared to RK6. The parameter list of RK6A 
differs from the RK6 parameter list. Explanations are given in the appendix. 
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#13: bahn4 :z (RK6A(magnetdipol (q , 1), q, [0, 1.82145, -0.595132, 0.0132261, 

0.0150255, 0.0746424, 0.0648287], 500, 0.2, [2, 3, 4], -8))u[2, 3, 4] 


#14: (RK6A(magnetdipol (q , 1), q, [0, 

0.0746424, 0.0648287], 10, 0.2 



1.82145, -0.595132, 0.0132261, 0.0150255 
, [2, 3, 4], -8)) U [5, 6, 7] 



#14 


Schon eine winzige Anderung einer Anfangsbedingung andert die Bahn nach kurzem Verlauf ganz- 
lich. bahn5 entwickelt fur vj.1 > 285 eine besonders interessantes Detail. Die Berechnung erfordert 
viel Zeit, denn die Schrittweite sinkt gegen Ende der Integration betrachtlich. In der Statusleiste zeigt 
RK6A den Rechenfortschritt an. 


A tiny change of the initial values changes the orbit completely. bahn5 shows an especially 
remarkable detail for vj.1 > 285. Its calculations needs much time because the step size 
becomes very small at the end of the integration process. You can follow the calculation pro- 
gress in the status bar (bottom left). 

#15: bahn 5 := (RK6A(magnetdipol (q , 1), q, [0, 1.82145, -0.595132, 0.0132261, 

0.0150255, 0.0746423, 0.0648287], 289, 0.2, [2, 3, 4], -10))u[2, 3, 4] 
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Nicht weniger reizvoll ist die Bahn eines geladenen Teilchens in einem oszillierenden elektrischen 
Quadrupolfeld. Ein solches Feld ist das Wirkungsprinzip einer „Paul-Falle“ oder „Quadrupole Ion 
trap“. Auch diese Anwendung liefert eine dreidimensionale Bahn. Lasst man den Kasten im 
3D-Fenster rotieren, so erkennt man die Eleganz dieser schonen Bahnkurve erst richtig. Ausdruck #16 
liefert die rechte Seite der Differentialgleichung, Ausdruck #17 berechnet die Bahn. Rechenzeit bei 
3 GHz CPU-Takt ist ca. 40 s. 

The orbit of a charged particle in an oscillating quadrupolar field is not less attractive. Such a 
field is the action principle of a "Paul-Trap" or "Quadrupole Ion Trap". This application deliv- 
ers a three dimensional orbit, too. Rotating the box in the 3D Plot Window gives an impres- 
sion of the elegance of this beautiful orbit. Expression #16 gives the right side of the differen- 
tial equation and expression #17 calculates the orbit. Calculation time is approximately 30s 
(3 GHz CPU-beat). 

ionfall(w, a, o, p) := APPEND(w , a-COS(w -o + p)-[w , w , - 2-w ]) 

#16: [5 7] 1 L 2 3 4j 

#17: bahn6 :i (RK6(ionfall(q, 1, 2.098823851, 0), q, [0, 0, 1, 2, 0.56, 0, 0], 0.1, 

4000))||[2 , 3, 4] 



4. Das Drei-Korper-Problem 

Das Drei-Korper-Problem liefert eine Vielfalt an Bahnkurven. Ich mochte drei von ihnen vorstellen. 
Die ersten beiden sind Losungen des eingeschrankten Drei-Korper-Problems, bei dem die Masse des 
dritten Korpers so klein ist, dass ihr Einfluss auf die beiden anderen vemachlassigt werden kann. Die 
sog. Trojaner im Gravitationsfeld von Sonne und Jupiter sind ein Beispiel dafur. Mit den Ausdriicken 
#18 und #19 erfolgt die Berechnung. Die Hauptmassen haben die Positionen [-0.9999,0] und 
[0.0001,0] (im Bild als rote Punkte eingetragen). 

Die Rechenzeit betragt bei 3 GHz CPU-Takt ca. 25 sSekunden. 
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4. The Three Body Problem 

The Three Body Problem delivers a variety of orbits. I'd like to present three of them. The 
first both of them are solutions of the restricted Three Body Problem: in this case the mass of 
the third body is so little that its influence on both other bodies can be neglected. The so 
called Trojans in the gravity field of sun and Jupiter are an example for this. Expressions #18 
and #19 do the calculation. The main masses have positions [-0.9999,0] and [0.0001,0] 

(- which are the red points in the graph). Calculation time is approximately 25 s. 

drei koerp(w, ml, m2, dl, d2, om, k 1 , k2) :i 
Prog 

#18: kl := - ml/C Cwj.2 - dl) A 2 + W|3 A 2) A (3/2) 

k2 :i - m2 / ((wj.2 - d2) A 2 + vu3 A 2) A (3/2) 

RETURN [wj4, w;5, kl-(w;2 - dl) + k2-(wj2 - d2) + om- (om-wj2 + 2-wj5), 


kl*W|3 + k2-W|3 + om- (om-W|3 - 2-W|4)] 


#19: 


bahn7 := (RK6(drei koerpCq , 0.9999, 0.0001, 0.0001, -0.9999, 1), q, [0, -0.6, 
1.03923, 0.33773, 0.195], 0.16, 1651))u[2, 3] 



Tm Zusammenhang mit dem Mondlandeprogramm der NASA hat Richard Arenstorf periodische Bah- 
nen fur Raumfahrzeuge im Erde-Mond-System untersucht {New periodic solutions of the plane three- 
body-problem corresponding to elliptic motion in the lunar theory, J. Different. Equ. 4 (1968), 
202-256). Auch ohne Kenntnis dieser Arbeit kann man solche Bahnen mit etwas Geschick selber fin- 
den. Hier ist ein Beispiel einer solchen Bahn, wobei die Erde-Mond-Bewegung zu einer Rotation mit 
konstanter Winkelgeschwindigkeit und konstantem Erde-Mond-Abstand vereinfacht ist. Die Erde ist 
bei [-0.0123, 0] zu finden, der Mond bei [0.9877, 0]. Man beachte, dass das Raumfahrzeug in der klei- 
nen rechten Schleife dem Mond sehr nahe kommt und sehr stark beschleunigt wird. 

(Rechenzeit zum Auswerten des Ausdrucks #20 ca. 15 s bei 3 GHz CPU-Takt.) 
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Richard Arenstorf investigated in connection with the moon landing program of NASA peri- 
odic orbits for space crafts in the Earth-Moon-system (New periodic solutions of the plane 
three-body-problem corresponding to elliptic motion in the lunar theory, J. Different. Equ. 4 
(1968), 202-256). Even without knowing this paper but being a little bit skilful one can find 
such orbits by himself. Here is an example of such an orbit - the earth-moon-motion is sim- 
plified to a rotation with constant angular velocity and constant distance between earth and 
moon. You can find the earth at [-0.0123, 0], and the moon at [0.9877, 0]. Please note that 
the space craft is approaching moon very close in the right loop and that is accelerated very 
much. (Calculation time for evaluating expression #20 is approx. 15s.) 


# 20 : 


bahn8 := (RK6(drei koerpCq , 
1), q, [0, 1.002, 0, 0, 


0.987722529, 0.012277471, -0.012277471, 0.987722529, 
-1.35661886266766633], 0.01, 1120))w[2, 3] 




(RK6A(d re i koerpCq , 0.937722529, 0.012277471, -0.012277471, 0.937722529, 1), q, [0, 1.002, 0, 0, 
-1.35661336266766633], 11.2, 0.001, [2, 3], -10)) u [4, 5] 

Ich lade Sie ein, das Geschwindigkeitsdiagramm (rechtes Bild) mit den Funktionen RK, RK6 und RK6A 
zu zeichnen und die Ergebnisse zu vergleichen. 

I’d like to invite you comparing RK, RK6 and RK6A plotting the velocity diagram (right figure). 

Die Umrundung des Mondes muss sehr kleinschrittig erfolgen, da sonst erhebliche Diskretisierungs- 
fehler entstiinden. Hingegen vertragen die groBen Schleifen in der linken Halfte des Koordinatensys- 
tems etwa 20-fach groBere Integrationsschritte bei gleichem lokalem Diskretisierungsfehler pro 
Schritt. RK6A ist darum fur dieses Anfangswertproblem viel besser geeignet als RK6. Der folgende 
Ausdruck #2 1 berechnet die Schleifenbahn noch einmal mit RK6A. Die Rechenzeit ist kiirzer, das Er- 
gebnis wesentlich genauer. 

Looping the moon must be performed in very small steps, because otherwise we would find 
serious discretisation errors. On the other hand the big loops in the left half of the system of 
coordinates tolerate 20-fold larger integration steps having the same local discretisation error 
per step. So RK6A is much more appropriate for this IVP than RK6. The following expression 
#21 calculates the orbit of the loop once more but using RK6A. Calculation time is less and 
the result is significantly more accurate. 
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#21: bahn8a := (RK6A(drei koerptq , 0.987722529, 0.012277471, -0.012277471, 

0.987722529, 1), q, [0, 1.002, 0, 0, -1.35661886266766633], 11.2, 0.001, [2, 
3], -10» U [2, 3] 

Aussagekraftiger als die obige Darstellung im rotierenden Bezugssystem ist die folgende Darstellung 
der Bahn in einem Inertialsystem. Dazu miisscn die Positionsdaten um den Ursprung gedreht werden, 
wobei der Drehwinkel proportional zurzeit ist. Dies leistet die Funktion d reh (a , Q) . Weil bahn8a 
aus den errechneten RK6A-Ergebnissen die Zeitinformation ausblendet, berechnet Ausdruck #23 die 
Bahn noch einmal, schlieBt aber die 1 . Spalte der RK6A-Ergebnismatrix mit ein. Wenn man Ausdruck 
#24 approximiert, erhalt man die Bahndaten im Inertialsystem. Das folgende Bild zeigt die Bahn. 

More meaningful than the presentation in the rotating reference system (as given above) is 
the following presentation of the orbit in an inertial system. For this purpose we have to ro- 
tate the position data around the origin with the rotation angle being proportional to time. 
Function dreh(a,Q) does the job. #23 calculates the orbit once more - including the first 
column of the RK6A-result matrix - because bahn8a from #21 does not show the time infor- 
mation. Approximating #24 gives the orbit data for the inertial system. The following figure 
shows the orbit. 
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#22: dreh(a, Q) := VECTOR 




•[a , a 1, i, DIM(a) 



- SIN(0-a ) 

C0S(0-a ) 
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J 


#23: balinSb := (RK6A(drei koerptq , 0.987722529, 0.012277471, -0.012277471, 

0.987722529, 1), q, [0, 1.002, 0, 0, -1.35661886266766633], 23, 0.001, [2, 
3], -10)) U [1, 2, 3] 


#24: di’eh(bahn8b , -1) 
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Man verfolge das Raumfahrzeug: Bei [1.002, 0] nahe dem Mond startet es, lauft gegen den Uhrzeiger- 
sinn auf einer Kepler-Ellipse um den Schwerpunkt Erde-Mond gut zweidreiviertel Umdrehungen weit. 
Inzwischen hat sich der Mond gut eindreiviertel Umdrehungen im gleichen Drehsinn um den System- 
schwerpunkt gedreht, Raumfahrzeug und Mond begegnen einander wieder. Wahrend des Vorbeiflugs 
am Mond wird der Geschwindigkeitsvektor des Raumfahrzeugs in kurzer Zeit um ca. 270° im Uhrzei- 
gersinn gedreht. Das Raumfahrzeug durchlauft eine kleine Schleife. AnschlieBend lauft das Raum- 
fahrzeug auf einer neuen Ellipse um den Systemschwerpunkt bis zur nachsten Begegnung mit dem 
Mond. 

We will follow the space craft: It starts at [1 .002, 0] which is close to the moon, and then it 
runs counter clockwise on a Kepler ellipse around the earth-moon centre of gravity about two 
and three quarters of rotation. In the meanwhile moon has rotated about one and three quar- 
ters in the same direction of rotation around the centre of gravity of the system. Moon and 
space craft are meeting again. During passing moon the velocity vector of the space craft is 
rotated clockwise about approximately 270° within a short while. The space craft makes a 
small loop. And then it runs on a new ellipse around the system centre of gravity until its next 
meeting with the moon. 

Der Sonderfall dreier Korper mit gleicher Masse soil das abschlieBende Beispiel sein. Bei geeigneten 
Anfangsbedingungen laufen die drei Korper auf einer gemeinsamen Balm, die einer Lemniskate ah- 
nelt, einander hinterher. Notwendig dafur ist: Gesamtdrehimpuls L = 0, Gesamtimpulsanderung p'= 0, 
folglich kann man p = 0 wahlen und auch Schwerpunkt s = [0, 0], Legt man die Startorter r x (0) = 

= [-1, 0], r 2 (0) = [0, 0] und r 3 (0) = [1,0] fest, hat man nur noch zwei freie Parameter, um die Bahn zu 
beeinflussen. Ich habe den Betrag u und die Richtung a der Geschwindigkeit der Masse im Ursprung 
als Parameter gewahlt. Nach einigem Probieren war die Bahn gefunden: u - 1.271 und a = 0.994. Es 
macht nichts, dass diese Werte noch nicht sonderlich genau sind. Im Gegensatz zum vorigen Beispiel 
ist die Bahn auch unter kleinen Storungen stabil. Ausdruck #25 liefert die wechselseitigen Anzie- 
hungskrafte zwischen den Korpem, Ausdruck #26 berechnet die Bahnen der drei Korper. 

Die Funktion pu3 dient dazu, dass man diese auf einmal zeichnen kann. Man benutze dafur den Aus- 
druck #28. Rechenzeit bei 3 GHz CPU-Takt ca. 10 s. 

The final example is the special case of three bodies of the same masses. Applying appro- 
priate initial values the three bodies follow each other on a common orbit which is similar to a 
leminiscate. Necessary for this is a total angular momentum L = 0 and a total change of the 
momentum p' = 0. So one can choose p = 0 and the centre of gravity s = [0, 0]. Fixing the 
start positions r_\ (0) = [-1 , 0], r 2 (0) = [0, 0] and r 3 ( 0) = [1 , 0], only two free parameters are 
remaining in order to influence the orbit: I chose the absolute value u and the direction a of 
the velocity of the mass which is located in the origin. After some tries I found appropriate 
values: u = 1 .271 and a = 0.994. It does not matter that these values are not very accurate. 

In contrary to the example above this orbit remains stable even under small disturbances. 
Expression #25 gives the mutual forces of attraction between the three bodies: expression 
#26 calculates their orbits. 


Function pu3 serves to plot them simultaneously. Please use expression #28. Calculation 
time is approximately 10 s. 
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dreigleich(w, klx , kly , k2x , k2y , kBx , kBy , dl , d2 , dB) := 
Prog 

dl := - ((w|7 - w|ll) A 2 + (w|6 - w|10) A 2) A (- B/2) 
d2 := - ((will - w|B) A 2 + (w|10 - w|2) A 2) A (- B/2) 
dB := - ((w|3 - w|7) A 2 + (w|2 - w|6) A 2) A (- B/2) 


klx :: 

= dl ■ (w^6 - 

W|10) 

kly :: 

= dl ■ (yji 1 - 

W|ll) 

k2x :: 

- d2 ■ (w^lO 

- w|2) 

k2y :: 

= d2 ■ (will 

- w|B) 

kBx :: 

= dB ■ (w|2 - 

w|6) 

kBy :: 

= dB ■ (w|B - 

w|7) 


RETURN [w|4, w|5, kBx - k2x , kBy - k2y, w|S, w|9, klx - kBx , kly - kBy, w|12, W|1B, k2x - klx , k2y - kly] 


#26: bahn9(u, «) := RK6 dreigleich(q) , q, 


0 , - 1 , 0 , 


1 

u -COS(cx) , 

2 


1 

u • SIN (a) p 

2 


1 1 

0, 0, u • COS (<x) , u.SIN(a), 1, 0, u • COS(cx) , u.SIN(<x) 

2 2 


\ 

, 0.04, 158 


pu3(v, a, m) := 

Prog 

#27: m := bahn9(v, a) 

RETURN [m j, x [2 p 3], m U [6, 7], mu[10, 11]] 

#28: pu3(l. 271, 0.994) 



Zum Schluss mochte ich noch die Anregung aus dem DNL 62, S.21ff, aufgreifen und auf das letzte 
Beispiel anwenden. Im DNL wird beschrieben, wie man mit einem Schieber aus einem Vektor von 
Punktkoordinatenpaaren ein Koordinatenpaar dergestalt selektiert, dass DERIVE veranlasst wird, den 
entsprechenden Punkt in einem 2D-Diagramm zu aktualisieren. Das Vorgehen fiihrt zu dem Eindruck 
als bewege sich der Punkt, wenn der Schieber bewegt wird. Um die drei oben genannten. Korper zu 
animieren, zeichnet man den Ausdruck #30. Zuvor muss im 2D-Graphikfenster ein Schieber fur den 
Parameter t mit Bereich 1 < t < 159, Schrittzahl 158, eingefugt sein. Wichtige Einstellung: Punkte 
nicht verbinden. Vor dem Hintergrund einer bereits gezeichneten Bahn sollte man eine andere Punkt- 
farbe und groBe Punkte wahlen. Dieser „Trick t4 lasst sich natiirlich auch bei alien vorigen Beispielen 
einsetzen. 

Finally I’d like to come back to the suggestion given in DNL 62, pp 21, and apply it to the last 
example. There is described how to use a slider for selecting a pair of coordinates from a 
vector of such pairs that DERIVE is induced to move the respective point in a 2D-plot. This 
leads to the impression as if the point is moving when moving the slider. For animating the 
three bodies from above one has to plot expression #30 after having introduced a slider bar 
for parameter t with 1 <t< 159, 158 intervals with point size large and points not connected. 
It is recommended to plot the orbit and the points in different colours. Of course, this “trick” 
can be used for all examples above, too. 
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#29: 


mv ( a , t) := [VECTOR (a 

L l|t,i 


i, 2), VECTOR (a 

2 1 1 1 i 


i, 2), VECTOR(a , i , 2) 
3 , t , i 



Im Mathematik-Unterricht des Gymnasiums (im Bundesland Bayern, wo ich unterrichte) werden Run- 
ge-Kutta-Verfahren kaum Verwendung finden und zwar schon deshalb nicht, weil Differentialglei- 
chungen im Lehrplan nicht erscheinen. Anders ist es in Physik. Der neue Lehrplan der 10. Klasse sieht 
vor, konkrete eindimensionale Anfangswertprobleme durchzurechnen. Ein Beispiel ist die senkrecht 
beschleunigte Rakete, die bei konstanter Schubkraft zeitlinear Masse verliert. Zur Berechnung wird 
das Eulerverfahren verwendet, das tatsachlich ein Runge-Kutta-Verfahren darstellt, wenn auch das 
primitivste. Auf diese Vorbereitung aufbauend ist es denkbar, dass Schuler in einem nachfolgenden 
Seminar mit Hilfe eines Fertigbausteins, wie RK6 einer ist, sich weitere, auch kompliziertere Bewe- 
gungsaufgaben erschlieBen. 

Fur Lehrer konnen RK-Verfahren Hilfsmittel zum Erstellen spezifischer Simulationssoftware sein. Oft 
findet man nicht die Software, die der eigenen Methode entspricht. Dann ist ein CAS mit angeschlos- 
sener graphischer Darstellung der Rechenergebnisse nicht selten ein geeignetes Instrument. Mit RK6, 
und besonders RK6A, ist nun auch DERIVE in der Lage Ablaufe zu simulieren, denen anspruchsvolle 
Anfangswertprobleme zugrunde liegen. 


5. Final remark 

In gymnasium math teaching (in the federal state Bavaria where I am teaching) Runge-Kutta- 
methods are hardly used. The reason is easy to find: differential equations do not appear in 
the curriculum. Things are quite different in physics. The new curriculum for the 10th form 
provides calculating practical one dimensional initial value problems. One example is the 
vertical accelerated rocket with constant thrust losing mass linear wrt time. The Euler-method 
-which in fact is the most elementary Runge-Kutta-method - is used for calculation. Based 
on this preparation one can imagine that students in a following seminar supported by a 
ready made “Black Box“ - like RK6 - are able to make more complicated motion problems 
accessible for themselves. 

Teachers can use RK-methods as tools for creating specific simulation software. It occurs 
often that one cannot find the software which meets one’s own teaching method. Then a 
CAS together with integrated plotting facilities is in many cases an appropriate tool. With 
RK6, and especially with RK6A, DERIVE is now also able to simulate processes based on 
ambitious initial value problems. 
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6. Anhang 

RK6A(u , y , yO , xend , hO , sel , tol ) hat folgende Parameter: u, y, yO wie bei RK. xend ist das Ende 
der Integration. Integriert wird von yOj.1 bis xend (und ein bisschen dariiber hinaus). hO ist die 
Schrittweite des 1 . Schritts. Der Wert ist unkritisch. Die eingebaute Schrittweitensteuerung sorgt fur 
die Anpassung an richtige Werte. sel ist ein Vektor, dessen Elemente Indizes auf den Vcktor y sind. 
Ausgehend von den indizierten Elementen von y wird der lokale Diskretisierungsfehler abgeschatzt 
und aus diesem die Weite des folgenden RK-Schritts abgeleitet. Ist der Funktionswcrt y eindimensio- 
nal, setzt man sel = [2], bei zweidimensionalem y bedeutet sel = [2,3], dass die euklidische Norm aus 
beiden Koordinaten verwendet wird. tol ist ein MaB fur die Genauigkeit der Berechnung. tol ist unge- 
fahr loglO des relativen Fehlers der Fortschreitungsrichtung. Sinnvolle Werte liegen zwischen -6 und - 
16. Noch kleinere Werte sind nicht sinnvoll, es machen sich dann Rundungsfehler des Algorithmus 
bemerkbar. Bei groBeren Werten kommt der Algorithmus (abhangig vom Anfangswertproblem) ins 
„Stottem“, es mussen viele Einzelschritte zuriickgenommen und mit kleinerer Schrittweite wiederholt 
werden. Empfohlen wird tol = -10 fur erste Untersuchungen eines Anfangswertproblems. Die Variable 
hprot dient dazu, die Schrittweitensteuerung zu uberpriifen. hprot ist ein Protokoll der Korrekturfakto- 
ren der Schrittweitensteuerung. Die Werte sollten nahe bei 1 liegen. Werte kleiner als 0,8 bedeuten, 
dass Schritte verworfen werden mussten. Kommt das haufig vor, sollte tol kleiner gewahlt werden. 

6. Appendix 

RK6A(u , y , yO , xend , hO , sel , tol ) uses the following parameters: u, y, yO have the same 
meaning as in RK. xend is the end of integration. Integration is running from y0j,l to xend 
(and a little bit further). hO is the step size of the first step. Its value is not critical. The built in 
control of the step size takes care for adaptation on appropriate values, sel is a vector which 
elements are indices on vector y. Starting with the indicated elements of y the local discreti- 
sation error is estimated and this is the base for setting the following RK-step. In case of a 
one dimensional function value set sel = [2], in case of a two dimensional y then set sel = 
[2,3] which has the meaning that the Euclidean norm of both coordinates is used, tol is a 
measure for the accuracy of the calculation which is approximately log 10 of the relative error 
in the direction of progress. Values between -6 and -16 make sense. Smaller values don’t 
make sense because of rounding errors of the algorithm. Taking greater values make the 
algorithm “stuttering" dependent on the given IVP. Many single steps must be taken back 
and then repeated applying smaller steps. I recommend tol = -10 for the first investigations of 
an IVP. Variable hprot serves testing the step size control. It is a report of the correcting fac- 
tors of the step size control. Its values should lie close to 1 . Values less 0.8 say that the 
steps had to be rejected. If this occurs often then you should choose a smaller value for tol. 

Comment of the editor 

I left Heinrich Ludwig's contribution in its original German version to the convenience of the 
many German speaking members of the DUG. For our English speaking friends I translated 
it into English. So don’t blame Heinrich for any mistakes in the English version. 

Josef 
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Runge-Kutta - Unveiled 

Josef Bohm, Wiirmla, Austria 

When I read Heinrich Ludwig's great contribution on the Runge-Kutta-Method I found myself in- 
spired enough to inform (again after many years) about this standard method solving DEs and systems 
of DEs numerically - and how this method has been treated by DERIVE (and other tools) so far. This 
is a short summary of the procedure (which can be found in many textbooks): 

d u x „ , . 

— = f x (t,u x ,u 2 ,...,u m ) 
at 

du ? r . 

—- = f2(t,u l ,u 2 ,...,u m ) 


, ^1 ’ ^2 ? ***? ) 

at 

for a<t<b with initial values u x (< a ) = a x , u 2 {a) = a 2 , u m {a) = a m . 


The classic RK-method of order 4 is given by 

y' = f(t,y\ a<t<b, y(a) = a 


K =hf(t i ,w i ) 

i r r( h AU 

k 2 =hf t i +-,w i + — 

V z z J 

1 7 r( h k 2 ^l 

k 3 = h f 

v z z ) 

k A =h f(t i +h + w. + k 3 ) 


then w i+x = w. + —{k x + 2 k 2 + 2 k 3 +k 4 ) 
6 

for i = 0 , 1 ,..., N -l with h = - — — . 


Heinrich Ludwig refers in his article to DERIVE s built in RK- function (contained in the Utility file 
ODEApproximation . mth, which replaced ODE_APPR . MTH from earlier DERIVE versions. 
DERIVE changed and improved from the DOS-versions via Dfw4 and DERIVE 5 to DERIVE 6 offer- 
ing now its great programming features - but the RK-routine did not change. It still consists of some 
auxiliary functions which are not easy to follow for students because they call themselves applying the 
ITERATES-command in the last step. 

It might be difficult for students and for not experienced DERIVE-users to discover the important 
steps of the RK-method more or less hidden in the DERIVE code. As you might know from other 
occasions the ITERATES-command often slows down calculation. So I had the idea to use the pro- 
gramming features to rewrite the “classic” RK-Method. 


See first the original RK-procedure how it is given in DERIVE's Utility file: 
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# 1 : 

# 2 : 

#3: 

#4: 

#5: 


[PrecisionDigits := 12, NotationDigits := S] 

cl + ( lim p) + 2-(c2 + c3) 
v-hj_ + c3 


RK_AUX3(p , 
RK_AUX2(p , 
RK_AUXl(p , 
RK_AUX0(p , 


v , u_, 
v, u_, 
v, u_, 
v, vO, 


cl, c2, c3) := 

6 

cl, c2) := RK_AUX3(p, v, u_, cl, c2, lim p) 

v-hj_ + c2/2 

cl) := RK_AUX2(p, v, u_, cl, lim p) 

v-hj_ + cl/2 

n) := ITERATES(u_ + RK_AUXl(p, v, u_, lim p) , u_, vO, 

v-hj_ 


n) 


#6: RK(r , v, vO, h, n) := RK_AUX0(h ■ APPEND([1] , r) , v, vO, n) 

#7: RK([x ■ (1 - y), 0.3-yO - 1)], [t, x, y], [0, 2, 1], 0.1, 10) 


And here is my program which does without any auxiliary function and without any ITERATES. 

rk4(r, v, vO, h, n, i, kl, k2, k3, k4, tbl) := 

Prog 

tbl := [vO] 

i := 1 

Loop 

If i > n 

RETURN tbl 

kl := VECTOR (SUB5T(k, v, vO) , k, r) 

k2 := 5UBST(r , v, ADJOIN(vO|l + h/2, REST(vO) + h/2-kl)) 

k3 := 5UBST(r , v, ADJOIN(vO|l + h/2, REST (vO) + h/2-k2)) 

k4 := 5UB5T(r , v, ADJOIN(vO|l + h, REST (vO) + h-k3)) 

vO := ADJOIN (vGil + h, REST(vO) + h/6 - [1 , 2, 2, 1] - [kl , k2, k3, k4 ] ) 

tbl := APPENDCtbl, [vO]) 

i :+ 1 

The parameter list for entering the function remains the same. 

I tested my program on several examples. The next one is an also “classic” prey-predator-system (one 
predator- and two prey populations). This example was provided by Josef Lechner (some years ago 
and based on a book on Systems Dynamic by Bossel, 1994). I will come back again to this model 
later. 

Compare the calculation times. (I select the last row of the result matrix to demonstrate that both func- 
tions are giving the same results.) 

(RK( [0 . 1* x - 0.1-x-z, 0.12-y - 0.1-y-z, 0.1-x-z + 0.1-y-z - 0.1-z], [t, x, y, 

#19: 

z], [0 P 1, 1, 1], 1 P 200)) 

201 

#20: [200 p 0. 010246976, 0.55946290, 1.9882390] 

needs 0.61 sec. 

(rk4([0.1-x - 0.1-x-z, 0.12-y - 0.1-y-z, 0.1-x-z + 0.1-y-z - 0.1-z], [t, x, y, 

#21: 


z], [0, 1, 1, 1], 1, 200)) 

201 

#22: [200, 0.010246976, 0.55946290, 1.9882390] 

needs 0.31 sec. 


In this case RK6 from Heinrich Ludwig does not give much better results - but calculation time in- 
creases up to 2.4 sec : [200, 0.010246757, 0.55945398, 1.9882254]. 
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I was satisfied with my program - but not satisfied enough to change the subject. I thought “And what 
about the TIs?” I browsed in my collection of papers and what did I find? A letter from our DUG- 
member from the very first days, Richard Schorn. In 1 996 he sent a diskette with a TI-92 program for 
treating the prey-predator-model using the “classic” Runge-Kutta-method. In his letter Richard re- 
minds on his comment in DNL#5 where he gave a numerical DERIVE solution applying the RK- 
function from the utility file. This is the header of Richard’s letter from July 1996: 

OSH) 1R, Richard Schorn 
Am B lei changer 10 
D-87600 Kaufbeuren 

Kaufbcurcn, 13. Juli 1996 

As you can see, (almost) nothing gets lost in my stack of papers. Sometimes I forget one or the other 
of my treasures provided by the DUG-members but then suddenly they appear and ask for being pub- 
lished. 

hafua ( ) 

Prgm 

Local kl ,k2,k3,k4,11 ,12,13,14 
0 . 04+a 
0.0015+b 
0.03+c 
0.0017+d 
1-+h 

Cl rGraph : Cl rDraw: FnOff 
'20^xmi n : 100^xmax : 10^xscl 
"20^ymi n : 100^ymax : 10^yscl : 2^xres 
setGraph ( "Axes" , "On") 
c/d^xs : a/b^ys 
PtOn xs,ys 
For k, 1 , 5 

0^t : xs+5^x0 : ys+7^y0 : ' 1 : x0^xs : y0^ys 
6^dd 

While dd>2 
i +1 -*-i 

abs (xs-x0)+abs (ys-y0)^dd 
If i <20 : 6^dd 
f (x0,y0)+k1 :g(x0,y0)-11 

f (x0+k1 /2,y0+11 /2)^k2:g(x0+k1 12 , y0+l 1 / 2)->l 2 
f (x0+k2/2,y0+12/2)^k3:g(x0+k2/2,y0+12/2)^13 
f ( x0+k3 , y0+l 3 ) -*k4 : g ( x0+k3 , y0+l 3)^1 4 
x0+ (kl +2*k2+2*k3+k4) /6^x0 
y0+ (1 1 +2*12+2*1 3+14) /6^y0 
t+h^t 

If mod (i , 3 ) =0 : PtOn x0,y0 
EndWhile 
EndFor 
EndPrgm 

Changing for plotting a-t-x or t -y diagram needs some change in the program code. The “ordinary” 
user will in most cases not be able to perform the appropriate changes. 

The problem given in DNL#5 was: 

RK([ (0 . 07-0 . 0051y)x , (-0 . 08+0 . 0012x)y] , [t , x , y] , [0 , 160 , 10] , 1 , 200) 

If you want to solve it using haf ua ( ) then you have to edit the program. This might be too compli- 
cated for some users. On the other hand, in 1996 the TI-92 was in its childhood and programming on 
the TI was really a challenge. Thanks Richard for facing this challenge and Hello from Wiirmla. 


fFiTTOT F ^ T Y F]:T Y F4t Y fe Y F6nr 
\ -r |fl 1 gebra |Ca 1 c |0ther |Prgm 1 0 |C 1 ean 


■ h (a x - b-; 

< y) * f(x ? y) 


Done 

■ h( -cy + d 

x y) ■* g(x ? y) 


Done 

■ haf uaO 



Done 


MAIN 

RAD AUTO 

FUNC 3/30 




The coefficients and the initial values 
are given within the program which 
produces a family of phase diagrams. 
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Now I felt another challenge: producing a more user friendly program. Program rk ( ) is the result of 
my efforts: 



ifiT^T FT? fht v— 

f— |Rlgebra|CalG|Qther|FrgmIQ|Clean 


ijciean UpT~ 1 


x 1 =f £t ? x ? y> : 
y l= g<t ? x ? y>: 

x 1 = £ . 07- . 0051 y>*x 

y'=< -,GS+.GG12x>+J 

, CEnter=OK } CESC=CflHCEL> J 


■r-kQ 

rkO 


Error: Break 


l^ 1 1 qebralcaYcToiherlprgn I pic 1 ealT Up 

n 

■ r 

■ r 

f 


initial oal par. : 
ini_oal 1st oar. : 
ini_oal 2nd oar. : 

G 

3ne 

?ak 

160 

1G 

step size: |l 

number of steps: |20Q 

£Enter=0K ) CESC=CFlNCEO 

rkO 

MAIN RAD AUTD FUNC 2/30 


The device is busy some time ... and 
then you are presented a pop up window 
to make your choice: All settings 

(WINDOW and PLOTS are set auto- 
matically by the program!) 



Press ENTER and the Pop Up Window will open again to receive your choice: 




It is quite better to have both time vs - graphs on the same axes because of different scaling. (You can 
interrupt the program by pressing the ON-key and then trace the plot.) But let’s have both populations 
in the same graph, and press the third option which shows both diagrams on the same axes: 
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DATA 








cl 

c2 

c3 

c4 

c5 


1 

0 

160 

10 




2 

1 

162.58 

11.204 




3 

2 

164. 12 

12.583 




4 

3 

164.43 

14. 149 




5 

4 

163.36 

15.902 




6 

5 

160.78 

17.834 




7 

6 

156.62 

19.919 





rlcl=0 


The last option gives the phase dia- 
gram. 

The values for t , x(t ) and y(t ) can be 
found in the data sheet rkdata. 


I don’t want to bore you by printing the program. You can find it among the files provided for 
downloading on our website. 

This was quite nice but I was not completely satisfied with myself. RK, RK4 and as well are much 
more flexible because they are not restricted to two equations only. Instead of improving my TI-92 / 
Voyage 200 program I downloaded the latest version of TI-Nspire (Version 2) and tried to write a 
program. 


The TINspire Challenge: 

I am starting with presenting some screen shots - treating one of Heinrich Ludwig’s examples with 
four equations (see page 11): 



rk() 

Please enter all necessary data between {}: 



What you will notice first is that something like a report containing all the input is presented. I’ll show 
how to do this. (It is a new feature of Nspire version 2). 
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rk() 


Please enter all necessary data between {}: 

Right sides of DEs in {}: {vx J vy,-x/(x A 2-y A 2) A (3/2),-y/(x A 2-y A 2) A (3/2)} 
Variables in{}, {par,...}: {t,x,y,vx,vy} 

Initial values in {}: {0,0. 2, 0,0, 3} 

Step size: 0.03141459 
Number of steps: 1 00| 


Values of the variables in same order as 
given in the variables" List are given in the 
rows of matrix tbl. The contents of its rows 
are given in lists 11, 12, ... for further use. 


I enter all data as requested be- 
tween {}. In this case it would be 
more comfortable to predefine the 
equations in the calculator eg 
eqs:={vx, vy, ..} because here you 
can use the tablets. In the dialogue 
boxes you cannot. (This should be 
improved!) 

The final result is a matrix tbl with 
the rows giving the approximated 
values of the variables appearing in 
the order of the 2 nd input. For plot- 
ting the diagrams it is necessary to 
have lists of values. The program 
does this job, too. It converts the 
rows of the matrix to lists. 


tbl 

o 

0.20000000 

0 

0 

3 


0.03141459 

0.18807955 

0.09237391 

-0.73515619 

2.82894771 


0.06282918 

0.15628744 

0.17543537 

-1.24496662 

2.44134419 


0.09424377 

0.11240027 

0.24552034 

-1.51593299 

2.02637756 


0.12565836 

0.06264693 

0.30338625 

-1.63276989 

1.66961718 


0.15707295 
0.01067245 
0.35116325 ► 
-1.66644438 
1.38317539 


II 

{0,0.G3141459,0.06282918,0.09424377,0. 12565836,0. 15707295,0. 18848754, 0.21990^ 
13 

{ 0,0. 09237391, 0.17543537, 0.24552034, 0.30338625, 0.351 16325,0. 39091079, 0.42427^ 
IS 

{3, 2. 82894771, 2. 44134419, 2. 02637756, 1.66961718, 1.38317539,1. 15607553, G.97465> 


Now you can produce the scatter plot and choose the attributes of your choice. 
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You may compare this figure with the DERIVE plot from page 12. 

Before showing the program I will present another nice example (provided by Josef Lechner a couple 
of years ago). Josef was and still is an expert in dynamic systems. The next screen shows the report of 
my input including text and dialog boxes as well. There are dialog-boxes with a text and a field to 
enter the requested input. The first line (“Please . . .”) is a “text box”. Both features are new! ! Explana- 
tions are given within the program code (page 32). The TI-Nspire program is running now. 


rk() 


Please enter all necessary data between {}: 

Right sides of DEs in {}: 

{x*(0.1 0000000-0.1 0000000*z),y*(0.1 2000000-0.1 0000000*z), 0.1 000000Q*x*z+0.1 0 
Q0000Qyz-G.1 QG0QGQ0*z} 

Variables in {}, {par,...}: {t,x,y,z} 

Initial values in {}: {0,1 ,1 ,1} 

Step size: 1 
Number of steps: 1 00| 




This is a nice coloured plot of the development of the three populations for the first 100 time steps. 
You can see three scatter plots on the same axes using lists 11 (time), 12 (prey 1), 13 (prey 2) and 14 
(predator). 
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What we cannot do with TI-Nspire but what we can do with DERIVE is producing a 3D-phase dia- 
gram of all populations! 



Believe it or not, the following is the graph of the x-y-z-phase diagram presented on the Nspire-CAS 
screen! I produced an isometric projection of the object onto the x-y-plane multiplying the matrix con- 
sisting of the three columns for x, y and z with the respective projection matrix (isometric projection), 
converted the resulting columns again into lists and plotted them as a scatter plot. 



p32 


Josef Bohm: Runge-Kutta - Unveiled 


D-N-L#77 


This is the TI-Nspire program. Comments are the lines starting with ©. There was one problem which 
made the program a lot bulkier than the DERIVE program: I miss my favourite DERIVE command, 
the powerful VECTOR-command which must be substituted by For-EndFor-loops. 

Define rk()= 

Prgm 

:Local kl,k2,k3,k4,I,j,rvO,df,dv 
:© The text and dialog boxes: 

:Text “Please enter all necessary data between {}:” 

:Request “Right sides of Des in {}:”,f 
:Request “Variables in {}, {par,...}:”,v 
:Request “Initial values in {}:”,v0 
:Request “Step size:”,h 
:Request “Number of steps:”, n 
:df:=dim(f):dv:=dim(v) 

: tbl : =newMat(n+ 1 ,dv) 

: tbl [ 1 ] : =list ► mat(vO) 

:i:=l 

:© The big loop for the n steps: 

:While i<n 

:rvO :=seq(vO [i] ,I,2,dv) 

:kl:=f: k2:=f:k3:=f:k4:=f 

:© The 4 loops for calculating the auxiliary values 
:© replacing the VECTOR-command in DERIVE 
:For j,l,dv 

:kl:=lim(kl,v[j],vO[j]) 

:EndFor 

:v2:=augment({v0[l]+((h)/(2))},rv0+((h)/(2))*kl) 

:For j,l,dv 

:k2:=lim(k2,v[j],v2[j]) 

:EndFor 

: v3:=augment( { vO [ 1 ] +((h)/(2)) } ,rv0+((h)/(2))*k2) 

:For j,l,dv 

:k3:=lim(k3,vQ],v3Q]) 

:EndFor 

: v4 : =augment( { vO [ 1 ] +h } , rvO +h*k3) 

:For j,l,dv 

:k4:=lim(k4,v [j ] ,v4 [j ] ) 

:EndFor 

:© Calculating the next values and filling the next row in the matrix 
:v0:=augment({v0[l]+h},rv0+((h)/(6))*(kl+2*k2+2*k3+k4)) 

:tbl [i+1 ] : =list ► mat(vO) 

:i:=i+l 

:EndWhile 

:tbl:=tbl T 

:Text “Values of the variables in same order as given in the variables' list are given in the rows of 
matrix tbl. The contents of its rows are given in lists 11, 12, ... for further use.” 

:© Using indirection to dynamically create the lists according to the number of variables 
:For I,l,dim(tbl)[l] 

:matHist(tbl[i]) -+ #(“l”&string(i)) 

:EndFor 

:EndPrgm 
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I can imagine that one does not intend to include programming into his mathematics course on one 
hand but that he wants to “unveil” the Runge-Kutta-Method as an improvement of the Euler Methods 
on the other hand because he/she is not satisfied with only using the Black Box RK/RK4 or RK6/RK6A or 
rk(). Spreadsheets or programs with spreadsheet features implemented like Nspire or GeoGebra are 
excellent tools for systems of two equations. In GeoGebra in its recent version it is not so comfortable 
because functions of two variables cannot be defined. I’d like to present the TI-Nspire realization 
combined with a nice application for the sliders to study the influence of the parameters and included 
with a moveable initial point. So finally we will have a very flexible model. 

This is a realization of the Prey-Predator-dynamic system. 

It is decribed by the system of differential equations: 

x' = ax-bxy = f(x,y) . . . x = number of prey animals 

y' = -cy+dxy = g(x,y) . . . y = number of predators 

The numeric solution is calculated using the Runge-Kutta-procedure with h = 0.1. 

All parameters a, b, c, d of the equations and the initial numbers of prey and 

predators can be controlled by sliders in the Graphs & Geometry Application. 

We show the phase diagram prey vs predators. 

The diagrams time vs pray and time vs predator could be added. 

In the Lists & Spreadsheet Application one can follow the 

Runge-Kutta-Approximation.| 


s 


* 

hQ 

1 

X 

Q 

II 

T 

Done 

g[x l y}.=-cyrdxy 

Done 

'll 

O 

0.10000000 


First of all we define the two equations in the Calculator. Then we switch to the spreadsheet: 



h 


h 


\ 


( h h \ ( 

F2:f{d\,e\\ H2:f dl + — f2,e\ + — g2 , J2:f dl + — h2,el + — i2 

V 2 2 J V 2 2 j 

In cells E2, G2, 12, K2 and M2 replace /by g. Then copy down row #2. 


L2:f(dl + h-j2,el + h-k2) 
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Cells cl and c2 are linked with the coordinates of the initial point. 


x' = a x - b x y 
y' = -c x + d xy 


a : = 1.6 


c := 1. 1 

— =r= 


2. 0. 


b : = .8 


d =16 x_start =10. y_start=5. 


.5 o 


P=rr 


= p= 
500 0. 



V- Ptpv 

1 

® a sl 

x-prey 

;y<-pred 



^ ;S2 

x«-x_s 

y<- y- S 


jjiK j 

m s3 

-i 

y<- 



x' = ax~bxy 
y = -c x + d xy 


a := 1. 


c := .7 
= —0= 


Pi'edator 



Prey 


d =.l 


x_start =10. y_stari =20. 


— rT 


500 0. 




The last realisation shows another form of animation: You can grab and drag the initial point and 
demonstrate the difference between Euler Method and Improved Euler Method. In a second plot win- 
dow you can show the time vs population diagrams. Grabbing the point could be implemented in the 
RK-model, too (instead of the sliders for its coordinates). 



References: J.F. Faires, R.L. Burden: Numerische Methoden, Spektrum, 1994 

Steven Schonefeld, Numerical Analysis via Derive , MathWare, 1994 
C.H. Edwards, D.E. Penney: Differential Equations , Prentice Hall, 1996 
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The Octahedron of Horror 

Wolfgang Alvermann, Emden, Germany 
(translated by J. Bohm) 

An octahedron is a regular polyhedron. Its surface is formed by eight congruent equilateral triangles. 
Each octahedron can be inscribed in a cube such that the vertices of the octahedron are lying in the 
centre points of the faces of the cube. 


The octahedron AB CDSjS 2 presented in figure 1 (Abbildung 1) is given by its vertices A( 13 | -5 | 3), 
5(1 1 | 3 | 1), C(5 | 3 | 7) and 5/(13 | 1 | 9). This octahedron is inscribed a cube with vertices Pj through 
P 8 as explained above. 



a) The distance between two parallel faces of the octahedron is called the “Thickness of the body”. 
Calculate the thickness of the given octahedron as distance between point C and plane ABSj. 

b) Find the coordinates of the vertices P 6 and P 8 of the cube given in figure 1 . 

c) The midpoint of segment AB is M AB , midpoint of segment CD is M C d • Let g the straight line con- 
necting the two midpoints. The octahedron is rotated around g as rotation axis in such a way that 
point ^4(13|-5|3) is moved to anew position A ’(1 2+2 V2 | -1+V2 | 2+2a/2). 

Show that the respective rotation angle is a = 90°. 

What are the coordinates of point 5’ as new position of point B after rotating the cube? 

d) A family of planes E a : 2xj + x 2 + 2x 3 + 9-(2a-5) = 0, cke R , is given; let h a straight line passing Sj 
and 5^(5 | -3 | 1). Show that each plane E a of the family is perpendicular wrt to line h. Find the in- 
tersection point P a of plane E a with line h. 

[Control: P a ( 13 - 4a \ 1 - 2a | 9 - 4 a)\ 

For 0 < a < 1 the plane E a cuts off a pyramid with vertex Sj from the octahedron (see figure 2). 
Find the volume V a of this cut off pyramid. 
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e) Six pyramids with equal volumes V a are cut off the octahedron in such a way that each vertex of 
the octahedron is the vertex of a pyramid and the base of each cut off pyramid is parallel to the 

f n 

opposite face of the cube (compare with task d)). A body R a 0 < a < — is remaining. 

V 2) 

Explain the properties of this remaining body R a for a = -j and a = ^ with respect to the number 
and properties of its side faces. 


It is recommended entering the given points 
as variables before starting calculation. 
ps2 will follow later. 

pa + (pc-pb) = pcf 

Definition of plane E AB si', calculation of the 
normal vector and the normal unit vector. 

Formula for the distance from the formulary 

d = \ 

AB = S,P 6 = -S,P 8 

So the edge length of the cube a = 12 
because of d = a j2. 

The edge length of the octahedron is 6- V2. 


Task b) 



Task a) 



Solutions 
Entering the data 
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Task c) 

Point M A b remains fixed during the rotation; it must be shown that the dot product of 

amZ-amZ= o . 



B’ is calculated from 

p^,=p;. + 2-(mZ-p;.) 

B' = (l2-2V2 1 1- V2 1 2-2V2) 


Task d) 

To show: The normal vector n 2 of plane E a and the direction vector v of straight line h 

are collinear; then follows that E 1 h. 



m = -4 confirms collinearity. 

Substitution of x-, y- and z-coordinates of the 
line into the equation of the plane leads to 




Substitution of / = — into the equation of the 
line confirms the given point. 

The height of the pyramid is 6a; calculation 
of its base edge is remaining. 









p38 


Wolfgang Alvermann: The Octahedron of Horror 


D-N-L#77 



Task e) 


a 


3 


6 squares and 8 regular hexagons 


The DERIVE sketch of the given problem 




Pictures of a polyamide - model 
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Comments: 

> Tasks a) and b) are standard problems in analytic geometry applying vector cal- 
culation. 

> Application of the dot product in task c) is difficult because it demands a high de- 
gree of spatial sense which cannot be expected from average students. 

> In my opinion task d) is too complex for an end examination problem (see Pro- 
fessor Steinberg's comment below). 

> Recognizing the remaining bodies makes problems not only for students! 

Comment of Prof. Steinberg (University Oldenburg): 

I treat the most beautiful and most difficult problems during my 
teaching, the less difficult I pose in tests and the easiest prob- 
lems are appropriate for end examinations (Abitur). 

I completely share his opinion, W.A. 

Reference: 

http://www.spieqel.de/schulspieqel/wissen/Q. 151 8.558559.00.html 

I could not resist modelling and plotting the two remaining bodies with DERIVE. The right 
figure was produced applying the POLYGON_FILL-function which is provided among the util- 
ity files. Josef 



About the author: 

studied Mechanical Engineering and Mathematics at TU Hannover 1970 - 1975, 

since 1977 teacher at Vocational Schools II Emden 

(Fachgymnasium Technik, Fachoberschule Technik, Fachschule Technik) 

Aim: Including real life applications in mathematics teaching 
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TIME 2010, MALAGA; Spain 

The Invited Lectures 

Barbel Barzel, Germany 

The role of CAS when learning algebra and developing functional thinking 

Eugenio Roanes-Lozano, Spain 

Mechanical Theorem Proving in Geometry with DERIVE using Wu’s Method 

Michel Beaudin, Canada 

Using the Real Power of Computer Algebra 

Colette Laborde, France 

Tangible and dynamic mathematics 

Lectures and Workshops (WS) accepted for the DERIVE 
and TI-Nspire strand so far 

C++ as a programming language for CAS 
Coding Theory for the Classroom (WS) 

20 Years International DERIVE and CAS-TI user Group 

Semiotic Acts and Process in the definition of Limit of a function according to Weierstrass 
Semiotics Tools and Mathematics Curriculum 

From Inductive Reasoning to Proof by Induction with Geometry Expressions (WS) 

Mathematical Modelling with Geometry Expressions and TI-Nspire CAS 
Teaching Linear Regression in the PC Lab 

Solving Second Order ODEs, Two Non-analytical Methods Revisited 
Using Rational Arithmetic to Develop a Proof, “What Josef and Carl Saw ” 

Beyond Newton’s Law of Cooling - Updated Methods for Estimating Time Since Death 
Using DERIVE®’s Graphics Tools for Geographic Profiling of Serial Offenders 
Using Mathematical Tools in Forensic Investigations (WS) 

Laplace Transforms, ODEs and CAS 

Could it be possible to replace DERIVE by MAXIMA? 

WIRIS - DERIVE, enabling compatibility 

Construction of mathematical knowledge using graphic calculators (TI-84 plus & CAS) 
in the mathematics classroom 

3D Analytic geometry in using the link between CAS and Geometry applications 
of the TI-NSpire (WS) 

WIRIS, a complete CAS with DGS (WS) 
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Experiments with geometric loci 
Mathematical Models in Biology (WS) 

The ElGamal Public-Key Cryptosystem - A Full Implementation Using DERIVE 
A virtual laboratory for blended-leaming: Numerical Methods using WIRIS 
Integrating a new DERIVE 6 Video User Guide into Virtual Teaching 
A CAS routine for obtaining eigenfunctions for Bryan’s effect 
Interactive Explorations of Mathematics with TI-Nspire Technology 
1 1 Years of Master Thesis in Engineering using DERIVE in the University of Malaga 
Teaching Differential Equations and its Applications Using DERIVE 6 as a PeCas 
On some images of a straight line mapped via the harmonic cross-ratio 
Teaching Calculus and Numerical Analysis using CAS according to Bologna Process 
E-Learning and Joomla 

Using Notes with Interactive Math Boxes in TI-Nspire Software Version 2 (WS) 

Exploring Rose Curves with the TI-Nspire Calculator (WS) 

Lectures and Workshops (WS) accepted for the ACDCA-strand so far: 

Teaching and Assessing Polygons Using Technology 

Learning Math, Doing Math: Fractals and Dynamic Geometry (WS) 

Enhancing Mathematical Problem Solving through the use of Dynamic Software Autograph 
Arguing and Reasoning supported by Technology 
Base Technologies for Tutoring 

The Positive Aspects of Modeling Process in Teaching Mathematics 
Sliders - A Dynamic Support for Teaching Mathematics 
Physics Through GeoGebra Window 

Web based education and assessment in the Bologna process 

Problems and Prospects of Remote Teacher Training in Uniform Environment of E-Learning 
Mathematics lessons and Classroom examples, inspired by the "Dresden Morning Post” 
Software GOLUCA: Knowledge representation in Mental Calculation 
WebQuest on Conic Sections as a Learning Tool for Undergraduate Students 
An analysis of arguments for and against the CAS 

Competence, didactic situations and Virtual Environments for Teaching and Learning 
The Intergeo Project (WS) 

Analogy and dynamic geometry software together in approaching 3d geometry 

Expanding Room for Tacit Knowledge in Mathematics Education 

Using Technology to Support Mathematics Teaching 

Teaching Math With Advanced Learning Blocks 

Flexible Mathematics Content Preparation 
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Designing Task for CAS Classrooms 

Web Based Multimedia Symbolic Package (Linear Motion) Using a Computer Algebra System (CAS) 

The use of notebooks in mathematics instruction. What is manageable? What should be avoided? 
A field report after 10 years of CAS-application 

The Role of Dynamic Geometry Software in the Process of Learning: 

GeoGebra Example about Triangles 

Visualization and Interface in Exploration of Functions 

The Development of Immunology Educational Package using Computer Algebra System 
(CAS) Modelling 

Categorizing CAS Use within One Reform-Oriented United States Mathematics Textbook 
Critical Issues of Technology Use in Undergraduate Mathematics 

On Assessment of Teaching A Mathematical Topic Using Neural Networks Models 
(with a case study) 

Utilization of Free Software on Teaching Integral Function in the Schools of Economics 
The Impact of Computer Use on the Teaching of Geometry in Grade 8 

Folding and unfolding cones and cylinders with Cabri 3D: how to do it and how 
to use it in the classroom 

3D-Dynamical Geometry in Building Construction 

Those Magic Moments when you realise you could not do this with chalk! (WS) 

Overcoming difficulties in understanding of the nonlinear programming concepts 
A Computational Measure of Heterogeneity on Mathematical Skills 
Using Maxima in the Mathematics Classroom 

Applications of Multimedia Technology to study the ordinal thinking evolution of scholars 
from 3 to 7 years old 

The Past and the Future of Computer Algebra and its Integration into Mathematics Education 
Technology and the Yin&Yang of Mathematics Education 

An example of learning based on competences: Use of Maxima in Linear Algebra for Engineers 

Learning Math in the context of European Space for Higher Education 

Lessons learned using Publishers’Web-based software to assess student work 

Using intelligent adaptive assessment models for teaching mathematics 

Taking advantage of Sherman's march 

GeoGebra Workshop for the Initial Teacher Training in Primary Education 
Problem solving using Geogebra 

Dynamic Applets for Differential Equations and Dynamical Systems 

There are some lectures more which are now involved in the reviewing process. You can find all full 
abstracts together with all Conference details at http://www.time20 1 0.uma.es . 


